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Neurogenesis has been the subject of active research in recent years and many authors have 
explored the phenomenology of the process, its regulation and its purported purpose. Recent 
developments in bioluminescent imaging (BLI) allow direct in vivo imaging of neurogenesis, 
and in order to interpret the experimental results, mathematical models are necessary. This 
study proposes such a mathematical model that describes adult mammalian neurogenesis 
occurring in the subventricular zone and the subsequent migration of cells through the rostral 
migratory stream to the olfactory bulb (OB). This model assumes that a single chemoattrac- 
tant is responsible for cell migration, secreted both by the OB and in an endocrine fashion by 
the cells involved in neurogenesis. The solutions to the system of partial differential equations 
are compared with the physiological rodent process, as previously documented in the litera- 
ture and quantified through the use of BLI, and a parameter space is described, the 
corresponding solution to which matches that of the rodent model. A sensitivity analysis 
shows that this parameter space is stable to perturbation and furthermore that the system 
as a whole is sloppy. A large number of parameter sets are stochastically generated, and it 
is found that parameter spaces corresponding to physiologically plausible solutions generally 
obey constraints similar to the conditions reported in vivo. This further corroborates the 
model and its underlying assumptions based on the current understanding of the investigated 
phenomenon. Concomitantly, this leaves room for further quantitative predictions pertinent 
to the design of future proposed experiments. 
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1. INTRODUCTION 

Neurogenesis has been the subject of active research in 
recent years and many authors have explored the phe- 
nomenology of the process, its regulation and its 
purported purpose in both the subventricular zone 
(SVZ) and the dentate gyrus of the hippocampus. In 
the adult SVZ, type-B astrocytes that are present 
close to endothelial cells act as slowly dividing neural 
stem cells, capable of generating a progeny of type- A 
neuroblast precursors [1] through a rapidly proliferating 
type-C cell intermediate. Type-A neuroblasts then 
migrate in a unique fashion to the olfactory bulb 
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(OB) along a path known as the rostral migratory 
stream (RMS), which is large (of the order of milli- 
metres) in comparison with the size of each individual 
cell (of the order of microns) [2] . These neuronal precur- 
sors migrate in chains along the RMS, which is 
ensheathed by astrocytic processes and outlined by 
blood vessels [3-5]. After migrating through the RMS, 
these type-A neuroblasts arrive at the centre of the 
OB and then move radially outwards. They further 
specify into either OB granule cells or OB periglomeru- 
lar cells, and become fully developed mature neurons 
that stain with NeuN but are otherwise indistinguish- 
able from OB cells [1]. Thus, to a first approximation, 
the stem cells differentiate linearly in a cascade 
throughout their complex spatial journey. 

During migration, type-A cells continue to divide 
and initiate neuronal maturation [6], but their rate of 
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Figure 1. (a) Illustration of the approximate physiological location of the migratory process, (b) Location of the computational 
domain with inset a typical BLI. (c) The three interconnected boxes that form the computational domain; the initial concen- 
trations of the chemoattractant factor /a and cell density are shown in grey within the olfactory bulb (OB) and 
subventricular zone (SVZ), respectively. The central black box indicates the region considered to be the centre of the OB, 
where in this simplifying model all type- A neuroblasts which arrive either undergo apoptosis or specify into adult neurons. 
RMS, rostral migratory stream. Scale bar shows unit length in the model, corresponding to 1.4 mm if dimensionalized. 



proliferation is drastically reduced and the cell cycle 
time is lengthened [7]. The migration is predominantly 
towards the OB, but many parts of the individual paths 
can point in other directions. In answer to the ques- 
tion of what guides migration, it is often theorized 
that differentially expressed migration factors might 
be responsible for steering migration, and a very large 
number of factors have been shown to influence the pro- 
cess: brain-derived neurotrophic factor [8]; growth 
factors such as endocrine growth factor [9], nerve 
growth factor [10] and fibroblast growth factor-2 [11]; 
polysialylated neural cell adhesion molecules, [12] and 
many other substances have been shown to have a 
mixed effect upon the entire neurogenic process [13]. 
It therefore appears probable that an attractant factor 
in the OB might direct migrating type-A neuroblasts 
towards the bulb. However, cells migrate successfully 
even when the OB is surgically separated from the rest 
of the brain [14], although the magnitude of migration 
is greatly reduced. Insertion of an explant culture of OB 
cells subsequently ameliorates the process, but insertion 
of any other piece of brain tissue does not, nor do many 
known chemoattractants [15]. This implies that attrac- 
tant factors are both present within the OB and 
secreted from the cells in either an endocrine or paracrine 
fashion, or are present within the wider environment; the 
exact mechanism remains currently unknown. 

Recent developments in bioluminescent imaging 
(BLI) have allowed direct in vivo imaging of neurogen- 
esis within rodents over a prolonged time period and the 
measurement of a number of tagged cells temporally 
throughout the process and spatially throughout the 
brain [8]. Properties such as cell migration speeds, pro- 
liferation, specification and apoptosis rates undoubtedly 
have a significant influence upon the results of such 
experiments and would prove troublesome to quantify 
in rodents. Furthermore, the specific mechanism and 
chemical factors responsible for the regulation of the 
process are at present largely unknown as outlined ear- 
lier. We therefore propose a simplified mathematical 
model, the results from which will, we hope, guide 
further experimentation. Theoretical and computational 
modelling of neurogenesis is not new. It has been used to 
study processes during development [16] as well as in the 



postnatal and adult brain [17,18]. Aimone et al. [16] 
review a suite of neurogenesis models, classifying them 
into 'abstract', anatomy independent models and 'bio- 
logical' models that encompass specific details of the 
anatomical location. The latter category has been used 
to study processes in the hippocampus [19-23] and the 
OB [17]. What these models have in common though is 
their focus on the functional effect of neurogenesis. 
In this study, we do not aim to consider all the details 
of the functional effect of the processes we model; instead 
the focus here is on providing a tool to interpret the 
experimental measurements. The agreement between 
the results of this model and the physiological process 
lends credence to the assumptions upon which the 
model relies, most notably that migratory neuronal che- 
moattractants are secreted both at a constant rate by 
both the environment and in an endocrine/paracrine 
manner by the migrating cells. 

2. MATERIAL AND METHODS 
2.1. Model and assumptions 

Spatially, the migration of cells through the RMS is 
a complex three-dimensional process [1,13,24], but to 
date quantifying BLI only has occurred in two-dimensional 
axial slices. We therefore consider that modelling the 
process in two dimensions at greatly increased computatio- 
nal and mathematical ease is an acceptable simplification. 
We assume that the computational domain can take the 
form of two interconnected boxes, the larger corresponding 
to the anatomically larger OB and the smaller to the much 
smaller SVZ, as illustrated in figure 1. Given the fact 
that we are modelling neurogenesis in the adult brain, 
domain growth can be neglected, which further simplifies 
the mathematical framework [25] . 

Figure 2 shows a schematic overview of the model 
presented in this study. It is generally acknowledged 
that type-B cells in the SVZ specify into type-C tran- 
sit-amplifying cells, which then become the migratory 
type-A neuroblasts. These type-A migratory neuro- 
blasts do not change their morphology further before 
either becoming adult neurons or undergoing apoptosis 
[13]. Although a small proportion (1.6%) of type-A cells 
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Figure 2. A schematic overview of the model and its terms, illustrating the linear differentiation cascade of the neuroblasts and 
associated terms, y and e are apoptosis rates; /3 are proliferation rates; a and £ are differentiation rates; 8 are diffusion constants; 
/and g are chemoattractants. 



seem to be able to become type-E endothelial cells [26] 
while in the SVZ itself, we do not model the density 
of the structural type-E cells but rather include this 
small loss of type-A cells within their apoptotic term. 
Because this is constant spatially within the model, 
this is tantamount to assuming that the rate of 
growth of supporting endothelial cells is constant 
spatially throughout the brain. Hence, we assume that 
every cell specifies in a linear way (B — > C — > A in the 
SVZ) and we denote the rate of specification by 
i= A, B, C with ni = r^(r, t); i= A, B, C being the 
cell density of type-A, B, C cells, r the two-dimensional 
position vector and t the non-dimensionalized time. We 
allow type-A, -B and -C cells to undergo mitosis and 
apoptosis at different rates, and we assume that these 
rates are not a function of position for each cell, apart 
from the type-A migratory neuroblasts. As type-A 
cells have been reported to proliferate much more 
slowly outside the SVZ than within it and a large pro- 
portion of the new cells arriving in the OB die [27], 
we take the type-A cells to have two different rates of 
mitosis with logistic growth constants (3 Ai within the 
SVZ, and fi Ao outside it as well as two different rates 
of apoptosis, (yA+ e ) n A within the central region of 
the OB and 7a elsewhere in the domain. 

We take type-B and type-C cells to be fixed in space 
and we let type-A cells move under the influence of a 
single chemoattractant factor, f A (r,t). This factor is 
assumed to be diffusible, taken up by the type-A cells 
at a rate kd/a^a? able to undergo some natural decay 
A/a? produced by the other species of cell present at a 
rate and produced by the OB environment at a rate 



a% exp< 



(x- xq) 2 - (y- y 0 ) 2 



ai 



where the location of the centre of the OB is (xo,y 0 ) and 
ai and a 3 are positive constants. As a significantly sim- 
plifying step, we ignore the complex radial development 
and migration of type-A migratory neuroblasts into gran- 
ule cells or periglomerular cells [28], and we assume that 
type-A cells either specify into adult neurons (with cell 
density n N ) or die [27] , once they reach a region surround- 
ing the centre of the OB. For computational ease, the 



central region is a square of side b and we define e(r) 
and £(r) as two related functions that are non-zero only 
within it; 6 is the increased rate of type-A cell apoptosis 
as described earlier, and £ = /me is the rate of mature 
neuron formation within the OB, with /jl being a positive 
constant. We further assume that the mature neurons are 
diffusible, have a slow rate of apoptosis (yNn N ; y^^C 
yi^N) and, as a model simplification, undergo taxis 
against a constant environmental factor g that is present 
in the form of a Gaussian to represent the radial 
migration found in the physical system. If we denote dif- 
fusion constants by 5^ and taxis coefficients by rj h this 
leads to the following set of equations: 
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d for \x — xq\ < b and \y — yo\ < 6, 
0 otherwise, 
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£(r) oc e (r) = /ie(r); 
Alto == { 



/3 Ai within the SVZ, 
/3^ 0 otherwise. 



(2.8) 
(2.9) 



The system of partial differential equations (PDEs) 
defined earlier was closed with appropriate boundary 
(zero-flux) and initial conditions. According to the 
BLI experiment, the initial concentrations of mature 
neurons, type-C and type- A cells were set at zero, 
with type-B cells and the factor /a having an initial 
two-dimensional Gaussian distribution in the centre of 
the SVZ and the OB, respectively. The model was run 
for 300 days according to the experimental set-up [8]. 

2.2. Implementation 

The variables in this system of differential equations 
represent concentrations, either of cells or of chemical 
factors. It is therefore vital that the numerical solutions 
to such a system are never negative, in particular 
because allowing negativity for any brief period of 
time may introduce unphysical oscillations and greatly 
affect the stability of the system. In addition, the 
method of solution used must ensure the conservation 
of mass and impose zero-flux conditions at the bound- 
ary of the computational domain to represent the 
constraints both of this model and of anatomy. The 
algorithm used to solve this taxis -diffusion -reaction 
system was that presented by Gerisch & Chaplain 
[29]. This implementation uses the method of lines 
(MOLs) to separate the spatial and temporal discre- 
tization of differentials over the computational grid of 
15 600 spatial points across the domain. The diffusion 
and reaction terms in equations (2.1) -(2.5) were discre- 
tized using the standard second-order central difference 
approximation and pointwise evaluation method, 
respectively, at the central point of each grid cell, 
both of which were found to be sufficient in terms of 
accuracy. As reported previously [29,30], the discretiza- 
tion of the taxis terms required the application of 
upwinding techniques with the non-linear van Leer 
limiter function to guarantee accurate, positive sol- 
utions for the system of ordinary differential equations 
(MOL-ODEs). For the time domain, the efficient 
numerical integration scheme ROWMAP was used, a 
ROW-code of order 4 with Krylov techniques for large 
stiff ODEs [31], which has an automated choice of 
step-size to ensure that the local error caused in each 
iteration remains below a user-prescribed tolerance (in 
this case 10~ 6 ), while keeping the computational cost 
as low as possible. 

2.3. Model parameters 

Parameter values were inferred (in non-dimensional 
form) from the BLI dataset. We stochastically gener- 
ated values for each of the parameters over a range in 
mag nitude (10" 6 -10 6 ) using a Mersenne twister algor- 
ithm [32] for the pseudo-random number generator, 
with the majority of the solutions investigated using 
generated values for each parameter that lie in the 
region (0,10). The reason for this choice of space was 
pragmatic: as the whole system is non-dimensional, 



scaling the magnitude of every parameter essentially 
redefines the coordinate axes. By (manual) comparison 
between solutions from different parameter sets (sto- 
chastic or otherwise), a guided search for physiological 
solutions was initiated, which ceased when a parameter 
space could be described that did not contain un- 
physical solutions [33]. Generating purely random 
parameters within the interval (0,1) resulted in physio- 
logical solutions twice in 180 runs (1%); guided 
stochastic searches in narrower domains increased this 
to 30-40%. Searches for new parameter values were 
ceased when 200 random samples of a space defined 
by a +10 per cent perturbation to a physiologically 
appealing parameter space (described later) resulted 
in entirely physiological solutions. Computational 
time for parameter sets leading to physiological sol- 
utions was around 1800 s on a single core of an Intel 
Xeon Harpertown 5500 processor, whereas solutions 
leading to steep variations, large gradients, pattern for- 
mation or quantities close to zero were computationally 
more expensive. Of the system of equations (2.1) -(2.5), 
it is apparent that equations (2.1) and (2.2) have 
algebraic steady-state solutions of the form 



n B (t 
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Pb ~ 7b - olc 



-a A + /3 C - y c 



(2.10) 



+ J (a A - Pc + 7c) 



4a c /3 c (/3 B - 7b - «c) 



(2.11) 



where the unphysical solution corresponding to n B = 0 is 
ignored. We therefore constrained the parameter space 
by requiring equations (2.10) and (2.11) to be greater 
than zero in order to avoid trivial solutions. To avoid sol- 
utions where neurons randomly diffused through the 
RMS rather than through undergoing taxis, we required 
rji > Other biological taxis models [30] frequently 
have the diffusion term three orders of magnitude smaller 
than the taxis one; hence this constraint is not particu- 
larly stringent in context. To check that these 
conditions were sufficient for physical solutions, a small 
number of spaces that did not satisfy the positive 
steady-state conditions but with larger initial conditions 
were explored, in order to allow the possibility for taxis 
within a long period of non-equilibrium behaviour; 
unphysical oscillations were observed as the only result. 

2.4- Parameter and model sensitivity analysis 

The sensitivity of the solution to small perturbations in 
the parameter set was investigated in a straightforward 
manner by considering the space defined by subjecting 
each parameter value to a +10 per cent perturbation, 
stochastically sampling this space 200 times and then 
solving the resulting systems. This showed that spaces 
yielding physiological solutions with plausible taxis 
times (i.e. non-zero and approximately within an 
order of magnitude of the experimental data) are not 
unique. As is common in biological models, one has to 
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consider the system as a whole and observe correlations 
and correspondence with experimental data; little can 
be gained by trying to accurately determine the 
values of a small range of individual parameters. This 
problem is well documented, and similar to Gutenkunst 
et al. [34], we defined a generalized parameter-space 
least-squares metric £ 2 that quantifies how the solutions 
change throughout parameter space relative to one 
particular set of parameter coordinates 0* as 

£ 2 ( 0 ) = i^f c CT^^' *> r) " Vsc{e ' *> r) ] 2(lr dt 

(2.12) 

Here t represents the time ordinate, 6 represents 
the parameter space coordinates, X b and Y b are the 
maximum extent of the computational domain in 
the x- and ^/-directions, T c the time up to which the 
simulation was run and A = 2 • T C N S is a normalization 
factor with N s being the number of species within the 
model and the sum extending over all species present. 
The value of £ 2 was computed for every physical (non- 
oscillatory) solution obtained via fast trapezoidal 
numerical integration (trapezoidal integration performed 
successively across each dimension). Because of the fact 
that solving the PDE system is computationally intensive 
and most traditional ways of creating a representation 
of £ 2 would require a large number of solutions, we 
elected to use a non-linear interpolation method to 
find approximations of £ 2 straightforwardly; an interp- 
olation method that does not require an evenly spaced 
grid of function value points was chosen to allow the 
calculation of £ 2 for all solutions obtained. In this way, 
the stochastic generation of parameter values represents 
a Monte-Carlo approach to finding £ 2 . The 'obtuse 
angle' interpolation method was used [35] and found 
to provide an adequate approximation to the highly 
non-linear £ 2 function. The validity of this approach 
was checked by comparison with results obtained using 
other methods of non-linear interpolation, such as 
the use of radial basis functions, which displayed some 
evidence of fitting bias. Furthermore, in order to investi- 
gate the error in the approximation to £ 2 , simulations 
were run with parameters randomly perturbed by + 10 
per cent compared with 0*, i.e. within the whole range 
that the interpolation function was trained on. These 
parameters were noted, and the differences between 
interpolated and computed metrics evaluated. This indi- 
cates an average relative error of approximately 50 per 
cent (m= 0.49, n= 10). Close to a data point, relative 
error is significantly better. 

As documented by Gutenkunst et al. [34], the Hes- 
sian matrix of all possible second-order derivatives of 
the function £ 2 , i.e. the matrix 

lJ 86,86, ' 

has eigenvectors that are the principal axes of the such 
ellipsoids and has eigenvalues that span many orders of 
magnitude. We obtained this matrix from the interp- 
olated approximation to £ 2 via standard numerical 
high-order finite difference approximations (as detailed 



in [36]) with the use of suitable small step sizes. The 
value of having this matrix as a function of 6 (rather 
than 0*) is that it could form a basis for a collective par- 
ameter fit in future studies [25]. 

2.5. Experiments 

The experimental set-up considered in this study is the 
same as that presented in the study by Reumers et al. 
[8] where BLI was used to track in vivo stem cell 
migration. In summary, neural progenitor cells were 
transduced with lentiviral vectors encoding firefly luci- 
ferase that was injected into the SVZ of mice. BLI 
was used to detect and quantify the progeny of these 
transduced cells in the OB over a period of 45 weeks 
post-injection. On day 7 post-injection, a detectable 
BLI signal was present at the site of the injection and 
from week 4, a signal from the OB could also be 
detected. A full report of the experimental set-up and 
the results is given in [8]. 



3. RESULTS 

Solutions to this system broadly fell into two categories: 
those that were physiologically plausible and those that 
were patently not. Those that we considered plausi- 
ble showed signs of a steady, controlled migration, the 
relatively quick equilibration of n B and nc to their 
steady-state values and a steady growth of mature neur- 
ons in the OB over time until a maximum value was 
reached. An example of such a set of solutions in com- 
parison to the mouse BLI [8] is shown in figure 3. 
Taking the non-dimensional unit of time to be days, 
there is a good agreement between the model predic- 
tions and the experimental in vivo results, in terms of 
both the qualitative behaviour of the cells and the 
timing of the process. 

As discussed in §2.2, a number of constraints were 
placed on the parameter space from which solutions 
were generated; we found that 8f& ^> 5 a?n , which is 
indicative of molecular as opposed to cellular diffusion. 
Furthermore, we found that the logistic growth con- 
stants for all three types of cells within the SVZ were 
usually within the same order of magnitude. Spaces 
with migration-related parameter values that were 
representative of those that might be expected from the 
literature — /3c ^> Pa,Pb — often resulted in solutions 
that were not unphysical, but where the taxis of type-A 
cells was small with respect to the effect of death 
terms. This either resulted in no neuronal maturation 
within the OB or in a quick transition to a steady-state 
constant flux of neurons through the OB, and not 
in the extended buildup over time that has been 
observed. Solutions deemed to be exhibiting unphysical 
behaviour showed either signs of exponentially growing 
oscillations or spontaneous symmetry breaking with 
the onset of divergent pattern formation. The electronic 
supplementary material, SI contains an overview of 
all the constraints on the parameter space needed 
to generate physiological solutions, along with the 
parameter set used as a reference in this study. 

The sensitivity to small perturbations in the solution 
illustrated in figure 3 and of others that were similarly 
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Figure 3. (a, b) The emergence of mature neurons in the system and the time evolution of the chemotactic factor / A f° r the set of 
parameter values /3 B = 7.567, /3 C = 9.441, /3 Az = 1.735, /3 Ao = 0.785, a c = 2.182, a A = 1.700, y B = 0.701, y c = 0.774, y A = 0.737, 
y N = 0.079, 5 A = 0.037, 8 N = 0.053, 5 /A = 1.955, 7? A = 0.134, t 1n = 9.136, k a = 6.045, k b = 1.276, k c = 2.292, k d = 9.284, A = 
0.422, ax = 1.000, a 2 = 0.500, a 4 = 1.000, a 5 = 0.461, 6 = 0.500, d= 1.921, ^ = 0.213. Initial conditions were n c , n A , 0; 
n B and / A initially obeyed a Gaussian distribution: f A in the OB with a maximum magnitude of unity and a standard deviation 
of half of the OB's side; n^in the SVZ with a maximum magnitude of 100 and standard deviation of 5/6 x half of the side of the 
SVZ. (c)(i) A comparison between the OB cell density of the model (n A + n^) and (ii) recorded in vivo bioluminescent photon 
fluxes of the OB reported in [8]. (ps -1 ); photons per second. 



close to the observed behaviour in vivo was investiga- 
ted straightforwardly by subjecting every parameter 
value to a random +10 per cent perturbation and 
solving the resulting system. For the solution illustra- 
ted in figure 3, 200 such perturbed systems were 
generated and solved, and we observed no unphysical 
behaviour in any of these results. While the equilibrium 
densities of cells predictably changed, the overall quali- 
tative properties of the solutions remained the same 
and even solutions with radically different equilibrium 
densities displayed similar behaviour, as illustrated 
in figure 4. 

In the sloppy parameter analysis, we observed multi- 
dimensional ellipsoids for the £ 2 function, as shown in 
figure 5 a. This quantifies how the results of the model 
change as a function of its parameters, and therefore 
allows one to understand which parameters have the 
greatest effect upon the model outcome. Directions in 
parameter space, where a small change in parameter 
values results in a large change in £ 2 , are referred to 
as 'stiff axes'; those, where £ 2 deviates little, are 
known as 'sloppy'. As illustrated in figure 56, the mag- 
nitude of the eigenvalues of the Hessian matrix showed 
a substantial variation by decades. The eigenvalues 
span some 10 orders of magnitude when normalized 
by the value of the largest, indicating that the sloppiest 
axes of the ellipsoids illustrated are more than 10 5 times 



as long as the stiffest. Furthermore, the values are 
spread out within this interval, thus indicating that 
there is no sharp distinction between 'important' and 
'unimportant' parameter combinations and confirming 
the sloppy nature of the model under investigation 
[34]. This sloppiness is manifest in our failure to 
change the qualitative behaviour of the solutions by 
adjusting individual terms (e.g. decreasing migration 
time by increasing the magnitude of taxis terms 77). 



4. DISCUSSION 

In this study, we propose a mathematical model describ- 
ing adult mammalian neurogenesis occurring in the SVZ 
and the subsequent migration of cells through the RMS 
to the OB. This model assumes that a single chemoattrac- 
tant is responsible for cell migration, secreted both by the 
OB and in an endocrine fashion by the cells involved in 
neurogenesis. The solutions to the system of PDEs were 
compared with the physiological rodent process as 
observed by BLI [8] and a parameter space was described 
for which the corresponding solution matched that of the 
rodent model. A sensitivity analysis was carried out 
showing that the described parameter space is stable 
under perturbation and furthermore that the system as 
a whole is sloppy. A large number of parameter sets 
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Figure 4. A direct comparison between the most and least neurogenic (a and b, respectively) solutions obtained through a + 10% 
perturbation to the parameter space illustrated in figure 3. While the cell density in the OB varies greatly in magnitude between 
the two spaces (upper and lower bounds on the shaded region in (c); unperturbed solution in black), the overall behaviour 
changes very little. 



were stochastically generated, and it was found that par- 
ameter spaces corresponding to physiologically plausible 
solutions generally obey constraints similar to conditions 
reported in vivo in the literature. 

We note that a steady decay from the maximum OB 
photon flux (which is proportional to the cell density) is 
observed in the BLIs, but not in our model; we propose 
that this is owing to cells within the rodent being 
tagged within the SVZ initially at one point in time, 
which will then eventually die [37]. This model tracks 
all neurons generated within the SVZ and therefore 
we expect to see no decline with time once a steady 
state has been reached. In addition, this model does 
not suffer from experimental problems such as stereo- 
tactic surgery scars, which have been seen around the 
site of the SVZ injection in the rodent images in 
figure 3 a. 



A number of simplifications have been made in the 
course of this study. A two-dimensional domain was 
used to model sub ventricular neurogenesis, which is a 
frequently performed simplification in many mathe- 
matical models of biological processes. Furthermore, 
given that the RMS is a predominantly planar system, 
we do no expect the results to fundamentally change 
when incorporating the third dimension. Another 
simplification lies in the numerical method used to cal- 
culate the parameter-space least-squares metric £ 2 , i.e. 
through the use of a non-linear interpolation method. 
The highly non-linear nature of the metric could lead 
to larger errors away from the interpolation points 
which might influence the exact value of the derived 
Hessian. The representation of the migration used 
in this study is also a simplification. It is difficult to 
find or derive appropriate constraints for migration 
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Figure 5. (a) Illustrative contours of constant model deviation as quantified by £ 2 (generalized parameter-space least-squares 
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slices are through the /3 B -/3c plane (logistic growth constants for B- and C-type cells) with all other values being the same as 
those in figure 3. Contour values are labelled; the centre of the figure corresponds to a minimum, (b) The normalized absolute 
values of the eigenvalues of the Hessian matrix in the same parameter space location. These eigenvalues span nearly 10 
orders of magnitude; a characteristic of sloppy systems [32]. 



parameters from the literature. The majority of stu- 
dies investigate the regulatory effect of a plethora of 
factors upon neurogenesis. The subsequently available 
constraints have either been derived from in vitro exper- 
iments, which can only study migratory behaviour on 
artificial surfaces, or from slices of the brain [38-40]. 
While the former might be far from the physiological 
environment, the latter could be incomplete when, for 
example, sources of chemoattraction are not present in 
the studied section of the brain. BLI now allows us to 
investigate this migratory behaviour in vivo for the 
first time. 

Owing to the complexity of the underlying biology 
(involving not only migration with chemoattractants 
either binding or not to the extracellular matrix, but 
also cell differentiation, proliferation and death) it is 
not straightforward to derive quantitative information 
without the help of mathematical modelling. This study 
is a first step in this direction as it tries to incorporate 
the underlying biological mechanisms that influence 
the experimental output. Future research will focus on 
the incorporation of different migratory mechanisms. 
The model proposed in this study contains many 'free' 
parameters for which no values were available in the 
literature. Directly measuring them would be difficult 
for some parameters and impossible for others, especially 
in an in vivo environment. Collectively fitting the 
parameters to experimental data often yields large par- 
ameter uncertainties, as is the case for our study. 
However, it has been shown extensively in the literature 
that collective fitting in biological models could yield 
well-constrained predictions even when it left individual 
parameters very poorly constrained [34,41,42]. The 



results presented in [34] suggest that sloppy sensitivity 
spectra, as observed by the substantial variation in eigen- 
value magnitudes in the model presented here, are 
universal in systems biology models. This prevalence of 
sloppiness highlights the power of collective fits and 
suggests that the focus of modelling efforts should be on 
predictions rather than on parameters. 



5. CONCLUSION 

We have shown how a relatively simple system of differ- 
ential equations can accurately model the biological 
process of subventricular neurogenesis and that the 
individual parameter values of this system are of 
equal importance in determining the properties of its 
solutions. We have found a particular parameter space 
(which is probably not to be unique) that is stable to 
perturbation and that yields results that match others 
in the literature. We believe that this model could be 
extended and reapplied, once more accurate infor- 
mation on the various parameters is available, in 
order to make quantitative predictions of neurogenic 
behaviour. Furthermore, the agreement between this 
model and the BLI data implies that at least one che- 
moattractant factor is secreted in an endocrine/ 
paracrine manner and is responsible for directing 
migration, as it is produced by all types of neural 
stem cells present. Whether this is physiologically the 
case remains to be confirmed experimentally. 
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